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1.  Introduction 

The  relationship  between  the  energy  of  an  explosive  source  and  the  amplitude  t)f  the  sei.smic 
waves  which  are  radiated  into  the  far  held  has  been  a  primtuy  interest  of  the  verihcation  progr;un  since 
its  beginning  (Latter  et  al.,  1959).  The  problem  is  made  difficult  by  the  fact  that  the  seismic  energy 
represents  only  a  small  fraction  of  the  toad  energy.  Most  of  the  energy  of  the  explosion  is  deposited 
within  the  elastic  radius  by  a  series  of  complicated  non-linear  proces.ses.  Given  that  the  wave  propaga¬ 
tion  problem  beyond  the  elastic  radius  is  essentially  solved,  the  primary  dilhculty  concerns  the  treat¬ 
ment  of  the  non-linear  region  surrounding  the  .source.  A  number  of  computer  ctxles  have  been 
developed  for  modeling  this  region,  but  they  are  fairly  complicated,  involving  hydrodynmnic  effects, 
shock  waves,  and  non-linear  equations  of  state.  Because  of  the  basic  numerical  approach  which  is  Idl- 
lowed  in  these  codes,  they  do  not  readily  provide  insight  into  questions  abtiul  which  p;a:imeters  are 
playing  critical  roles  in  determining  the  radiated  ela.stic  waves.  This  has  motivated  tlie  investigation  of 
an  alternative  method  of  modeling  this  region  immediately  surrounding  an  explosive  source. 

The  basic  objective  of  the  reseiuch  described  in  this  report  is  to  explain  the  energy  tliat  is  radi¬ 
ated  into  the  elastic  region  by  an  explosive  .source.  Tlie  energy  which  is  deposited  within  tlie  clastic 
radius  is  not  of  direct  interest  in  the  sen,se  that  it  does  not  proptigate  as  fttf  as  the  elastic  region.  How¬ 
ever,  it  is  necessary  to  obtain  an  estimtue  of  the  touil  tunouni  of  Uris  energy  which  is  deposited  within 
the  elastic  radius,  as  this  must  be  subtracted  from  the  energy  of  the  explosion  in  order  to  detcrinine  die 
energy  that  reaches  the  elastic  region. 

The  basic  approach  being  investigated  is  to  use  ;ur  equivalent  elastic  trctiunent  for  tlie  region 
between  the  original  cavity  radius  and  die  elastic  radius,  fliis  concept  of  an  equivalent  elastic  iiiedium 
has  been  used  quite  successfully  by  earthquake  engineers  to  model  the  non-lineiu'  behavior  ol  soils  iliat 
occurs  during  strong  groutid  motion.  An  adviuitago  of  this  approach  is  tliat  it  is  possible  to  nuike  u.se  of 
the  analytical  .solutions  which  lue  available  for  the  linear  problem,  fhe  central  idea  of  the  nietluKl  is  to 
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make  the  material  properties  a  function  of  the  stress  in  the  outward  propagating  pressure  pulse  ruid 
obtain  the  results  in  the  form  of  a  simple  numerical  propagatit)tt  of  antilytical  solutions.  ITie  present 
formulation  relates  density  and  bulk  elastic  properties  to  tlie  petik  pressure  in  the  pressure  pulse  and 
shear  and  anelastic  properties  to  the  maximum  shem  stress.  The  material  prtrperties  ;ue  adjusted  in  an 
iterative  process  so  that  appropriate  values  are  pre.sent  in  die  vicinity  of  the  propagating  pressure  pulse 
While  this  approach  is  only  an  approximation  to  calculations  with  hydrodynamic  equation-of-state 
codes,  it  has  the  advantage  of  providing  simple  analytic  results  in  which  the  role  of  vtirious  nuxlel 
parameters  is  easily  investigated.  This  is  imptirtant  when  one  wiuits  to  conduct  a  .sensitivity  analysis 
over  a  wide  range  of  explosion  sizes  and  material  ptuameters. 

There  are  a  number  of  lines  of  evidence  which  suggest  that  iliis  type  of  approach  could  possibly 
provide  some  u.seful  results.  The  study  of  Denny  and  Johnson  (1991)  .seemed  to  indicate  that  relatively 
simple  scaling  reiatioaships  could  explain  a  wide  variety  of  explosion  data.  A  .single  settling  relation¬ 
ship  was  found  to  be  reasonably  successful  in  reconciling  explosion  data  that  spanned  a  broad  riutge  of 
explosion  yield  and  source  medium.  The  data  set  ;dso  spanned  different  tyiu's  of  explosions,  including 
nuclear  explosions,  chemical  expUtsions  in  the  held,  and  .small  chemical  explosions  in  Uie  laboratory. 
These  resutis  could  be  interpreted  to  metui  tliat.  if  the  explosion  data  c:ui  be  exphiined  by  simple  rcUi- 
lion.ships,  then  a  corre.spondingly  simple  ireaunciu  of  the  physics  of  the  problem  might  also  be  possible. 
There  are  a  number  of  examples  in  continuum  mechanics  where  a  general  treaunent  of  a  problem  which 
lumps  together  a  number  physical  meclKuii.sms  is  po.vsible.  Anelastic  attenuation  of  seismic  waves  is 
one  such  example,  as  the  basic  physics  of  the  prextess  is  not  fully  understtHtd  ;uid  it  is  likely  that  a 
number  of  different  processes  are  re.sponsible,  and  yet  it  is  ptwsible  to  describe  all  t)f  these  eflecis  with 
a  single  parameter,  the  quality  factor  Q.  Finally,  if  Uiis  problem  is  viewed  in  terms  of  energy,  it 
becomes  a  ba.sic  trmsport  problem,  with  the  prim.-ay  Uisk  that  of  determining  bow  much  of  die  energy 
is  absorbed  and  how  much  is  transported, 

2.  Equivalent  Elastic  Method 

The  region  immediately  .surrounding  an  explosive  .source  experiences  sirc.s.scs  of  suflicieni  magni¬ 
tude  St)  that,  at  lea.st  lemponuily,  Uie  material  properties  are  significantly  ditlercnt  than  in  llic  normal 
state  of  low  slre.ss.  The  basic  idea  of  the  equivalent  clastic  inclJiod  is  to  retain  the  assumption  ol  .1 
UnCc'u  elastic  consiiiuiive  relationship,  but  to  ailjust  the  elastic  constants  aeet)riiing  to  the  level  ol  the 
stress,  'Hie  process  is  dilierent  tor  the  hulk  and  shear  pro|x.'rties.  It  shoultl  he  poinicti  out  iliat  die 
material  being  considered  liere  is  actually  a  linear  vi.sco-elastic  material,  as  the  consiiiiitive  eqiiaiions 
me  assumed  to  have  both  elastic  and  viscous  terms. 


Consider  first  the  bulk  properties.  Majtenyi  and  Foster  (1992)  have  suggested  tui;ilyiic;il  expres¬ 
sions  for  describing  the  changes  in  density  and  compressioiuii  veltKity  caused  by  the  passage  of  a  one¬ 
dimensional  pressure  wave.  These  have  been  mcxlilied  bek)w  to  make  them  more  appropriate  lor  die 
case  of  waves  in  three  dimensions.  Let  p„  be  the  ultimate  or  maximum  density  and  P  be  the  pressure 
Then  the  dependence  of  density  upon  pressure  is  expressed  as 

p(^)  =  Pu  -(Pu-P»)  e"'*'  (2.1) 

where 


1  P« 

A  s  - - 

Pu  Pf) 


(2.2) 


and  Po  and  Jl„  are  low  pressure  values  of  the  density  ;uid  bulk  modulus,  respectively.  In  terms  ol  die 
low  pressure  values  of  the  P  velocity  and  S  veliK’ity 


The  compressionat  vekx:ity  as  a  function  of  pressure  is  given  by 


VJP)  = 


Pu 


P» 
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(2.4) 


Next  consider  the  .shear  properties.  In  this  case  Uie  .shem  modulus  and  the  shear  quality  factor  are 
adju.sted  according  to  the  shear  sU'e.ss.  However,  it  is  convenient  to  express  the  relationships  in  tenns  of 
shear  strain  rather  than  shear  stress.  Let  T|„a,  be  the  shear  .strength,  the  maxiinu/n  stie.ir  stress  tliat  die 
material  will  sustain,  and  define  the  reference  strain  by 

Here  Po  is  stress  sheiu  modulus,  which  can  be  obtained  from  die  low  stress  slietir  veliK'ity  by 

=  PuVi  (2.6) 

The  appitfent  suain  associated  with  the  shear  stress  x  is  delined  to  be 
T 

Then  the  simple  hyperbolic  stfess-strain  relation.ship  (Hiudin  and  Dmevich.  I't72h)  yicULs  the  clleciive 
she.ar  modulus 


p(e)  =  ji,, -—J— 

I  +  e  /('r 
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Similarly,  the  effective  djunping  ratio  is 


d(e)  = 


e/er 


1  +  eler 

where  d^^  is  the  maximum  damping  ratio  for  the  material. 

A  modified  hyperbolic  relationship  can  be  obtained  by  defining  a  hyperbolic  strain  as 


=  r 


l+a  e  *' 


where  a  and  b  are  material  dependent  constants.  Then 
1 


^l(e)  = 


1  +  e/, 


and 


(2.8) 


(2.9) 


(2.10) 


d(e)  =  d„ 


eh 


1  +  eh 


(2.11) 


Note  that  the  constants  a  and  h  used  in  defining  e*  may  be  different  for  die  calculatitin  of  n  and  d. 


3.  Linear  VLscu*Elu.stic  Solid 

An  important  aspect  of  die  equivalent  elastic  treaunent  as  it  is  being  used  in  this  problem  is  die 
treatment  of  energy  absorption.  This  is  achieved  by  treating  the  material  as  a  linear  visco-elastic  solid. 
Because  the  energy  ab,sorpdon  in  this  problem  is  quite  significant,  it  is  worth  reviewing  the  basic 
dehnidons  of  the  quantides  involved.  The  stxuting  point  is  to  assume  that  the  Boliziniui  principle  of 
superposition  can  be  applied  tuid  that  the  material  can  be  de.scribed  by  a  linctu-  hcrediitu^  constitutive 
equadon.  Then  the  stress-sUain  relationship  c;m  be  written  as 
/ 

Tj^(t)  =  J  de^iiy)  t.Tl) 

where  C,^h(1)  are  now  die  genertdized  vi.sco-elastic  moduli.  In  die  frequency  domain  this  becomes 

T,,((o)  =  C,/*/((n)  et,(o))  (.T2) 

where  the  following  I'ourier  translorin  pair  has  been  introduced 

C„,,(co)  —  ^C,,„(0  (.T,^) 


4 


In  order  to  simplify  the  algebra,  consider  only  one  of  the  elastic  mtxluli  ;md  then  the  basic  heredi¬ 
tary  constitutive  equation  can  be  written  as 


x(0  =  JC(r-y)  i/elY) 


(3.4) 


In  general  C{q))  will  be  a  complex  quantity  so 
C(®)  =  Re{C({0)} +/■  Im{C((o)) 


=  lC(a))l  e 


1*0/' 


where 


tan(0c)  = 


Im{C(m)) 


(3.5) 


(3.6) 


Re(C{C))} 

This  study  follows  a  practice  common  in  seismology  and  introduces  the  quality  factor  Q  through  the 
equation  (Aki  and  Richards,  1980). 

C(o}>  =  Re{C(o)))(l  +  i  (?-')  (3.7) 

so  that  the  inverse  quality  factor  is  defined  as 


c-‘  . 

Re{C(w)} 


(3.8) 


A  critical  difference  between  elastic  and  vi.sco-ela.stic  materials  invrtlves  the  the  fact  ih;it  energy  is 
not  conserved  for  the  latter.  For  a  general  deformation,  the  change  in  internal  strain  energy  is 


p[e(r)  -  £{0)]  =  J  T{y)-j-e{y)  <iy 


0 

Consider  this  result  for  the  case  of  a  simple  hiumonic  stress  held 
t(r)  =  Re{-i  =  w-rinfor) 


(3.9) 


(3,10) 


where  is  a  scalar  constant  giving  the  mtixiinum  stress.  For  the  associated  strain  the  correspondence 
principle  gives 


e({)  =  Re(-r 


T(t) 


C((i)) 


(3.11) 


In  the  course  of  a  single  cycle  the  clastic  siniin  energy  will  incrctise  and  ilien  decrea.se.  It  can  be 
shown  that  the  muxiinuin  strain  energy  is  readied  when 

3ii 


(or  =  ^  +  0,; 


(3  12) 


5 


and  die  maximum  value  is 


max  [pel 


tm.. 


,  3it 


2  -Re{C((i))](I  +  (~^+0c)f^*^(6c)) 


(3.13) 


At  the  end  of  a  single  cycle  (f  =  2]t/co)  there  will  be  a  net  loss  of  energy  given  by 

A[pel  =  7t  Im{C((0)l  (3.14) 

A  quantity  of  interest  is  the  ratio  of  the  energy  loss  per  radian  to  the  maximum  energy,  which  is  often 
used  to  define  a  quality  factor  directly  related  to  energy. 

Afpe} 


Q 


-1 


2jt  max  [pel 

Using  the  results  derived  above,  this  becomes 
_ Im(C(io>} 


(3.1.3) 


Q 


-1 


Re{C((o))(I  +  (■y+ec)tan(0£:)) 


(3.16) 


(3.17) 


If  the  viscous  effects  are  small  compared  to  the  elastic  effects,  then  (2"‘  <  1  and 
Q-‘  =  Q-' 

and  the  difference  in  the  two  definitions  of  Q  is  not  iinportimt.  However,  when  viscous  eflects  iue 
appreciable  the  difference  between  the  two  defiiiitions  must  be  taken  into  account.  ;uid  the  exact  rela¬ 
tionship  is 

.  Q-' _ 


Ql 


I  +  ( Y+uur'(CJ''))(?"‘ 


(3.  IS) 


In  experimental  work  it  i.s  common  to  define  the  damping  ratio  as 


d  s 


Afpe] 


2jr  maxltfr)]  max[e(f'l 


(3.1‘J) 


and  this  leads  to 

d  =  -jsin(0c) 

In  terms  of  the  seismic  dclinilUm  of  ()"*,  die  relationship  for  the  dtunping  ratio  is 
2d 


(3.201 


Cr'  - 


Vl-w- 


(.3.21) 
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shear  stress.  K/PA 


shear  strain,  percent 


Figure  1.  A  hysteresis  curve  and  the  parameters  which  can  be  used  to  describe  it.  The  area  of  the 
ellipse  is  a  measure  of  A[pel  and  the  slope  of  a  line  through  the  end  point.s  of  the  curve  is 
Sh  =  tan(()>).  The  maximum  values  of  stress  and  suuin  are  also  shown. 
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A  convenient  method  of  analyzing  the  energy  relationships  of  visco-elastic  materials  is  in  terms 
of  hysteresis  curves.  A  cyclic  loading  of  the  material  results  in  a  stress-strain  curve  such  a.s  that  shown 
in  Figure  1.  For  a  linear  visco-elas:.c  material  this  curve  is  an  ellipse.  The  area  enclosed  by  tlie  curve 
is  a  measure  of  the  energy  los.'^  jier  cycle  A[pel  and  the  maximum  stress  and  maximum  strain  are  al.so 
easily  measured  from  the  curve,  so  equation  (3.19)  can  be  used  to  estimate  the  damping  ratio  of  the 
material.  Ther  equation  (3.21)  leads  to  an  estimate  of  the  seismic  quality  factor. 

4  *  2  slope  of  a  straight  line  through  the  end  points  of  the  hysteresis  curve  can  be  used  to  provide 
an  estimate  of  the  modulus  of  the  elastic  constant  iClto)!.  If  5*  is  the  slope  of  this  line,  then  it  c;m  be 
shown  that 


lC((i))l  = 


[(.S,-^l)--4.‘>„-Q-^l'^  ^  .sV-1 
25,[l+(2--l'" 


Figure  1  illustrates  how  an  estimate  of  the  ehusiic  mtxlulus  can  be  obitiined  from  a  hysteresis 
curve.  The  four  quantities  max(e],  max[T],  (Ji,  and  A[pe]  erm  be  measured  directly  from  the  curve. 
Then  equation  3.22  yields  an  estimate  of  the  absolute  value  of  the  elastic  mcxlulus  and  equatiims  3.19, 
3.21,  and  3.8  yield  an  estimate  of  its  pha.se.  Hence,  the  complete  complex  ela.stie  modulus  can  he 
obtained  from  the  hysteresis  curve.  This  estimate  of  the  elastic  modulus  will  in  general  be  a  function 
of  the  level  of  stress  and  strain. 

Hysteresis  curves  such  as  shown  in  Figure  1  provide  a  convenient  inetlitxl  for  the  e.xperimenial 
determination  of  the  elastic  properties  of  etirth  materials.  Both  stress  and  su;un  ;ue  measured  ;ls  a 
cyclic  stress  is  applied  and  the  re.sult  is  a  direct  determination  of  the  hysteresis  curve.  In  general  tire 
resulting  curve  for  actual  earth  materials  will  not  be  an  ellipse,  but  the  quantities  ma,x[el,  m;e((T].  (>. 
and  A[pE]  can  still  be  extracted  from  the  curve  and  convened  to  :m  equivalent  line;u  elastic  nuxluius 
(.see  for  example  Idriss  and  Seed  (1908)  ;u>d  lUudin  ;uid  Dmcvich  (1972a)). 

Figure  2  shows  hysteresis  lixips  that  were  calculated  for  iluce  dilfcrcnt  stress  levels  in  a  material 
having  the  properties  of  wet  luff.  The.se  properties  were  taken  to  be  p„  =  2..S8  gm/eny,  e,  -  10~*.  ;uid 
-  0.30.  The  equivtJcni  ela.siic  treauneni  of  .section  2  w;ls  u.sed  to  adjust  the  material  propenies  as 
the  peak  sire.ss  and  strain  were  cluuiged.  In  tJie  case  of  compre.ssiomd  stress,  there  is  a  slight  increase 
in  the  real  ptut  of  the  modulus  as  die  .stress  incretuses  ;uid  tlie  imagiieuy  p;irt  also  increases.  In  the  case 
of  .she;ir  .stre.s.s,  there  is  a  mtirked  decrease  in  the  real  pitft  of  the  intxlulus  ;md  :ui  increase  in  the  ima¬ 
ginary  part  as  the  level  of  the  stress  increa.se.s.  The  numerical  vidues  corresponding  to  these  hysieiesis 
curves  ;irc  listed  in  Table  1.  In  this  table  case  0  is  die  low  sness  situation,  ;uid  ca.ses  I  tlinnigli  .1  aie 
tlie  situations  shown  in  I'ieure  2  as  tlie  stress  level  is  increased. 


Table  I.  Equivalent  elastic  consumts  lor  liie  hysteresis  cur'es  of  Figure  2. 


Case  0 

Case  1 

Case  2 

C.ase  3 

Peak  pressure  (MPA) 

0 

5 

10 

30 

Peak  shear  stress  (MPA) 

0 

1 

3 

10 

Density  (gm/cc) 

2.10 

2.10 

2.10 

2.11 

P  velocity  (km/sec) 

2.85 

2.85 

2.85 

2.86 

S  velocity  (km/sec) 

1.45 

0.94 

0.69 

0.41 

e;‘ 

0.01 

0.18 

0.23 

0.26 

(?r‘ 

0.02 

0.48 

0.63 

0.72 

ea 


4.  Elastic  Wave  Propagator  Solution 

The  use  of  the  equivalent  elastic  method  requires  that  the  clasUKlyiuunic  equations  of  motion  be 
solved  for  an  explosive  source  in  an  inhomogeneous  medium.  Si:uid;u-d  propagaioi  mcihisjs  c;in  be 
used  to  solve  this  problem. 

Consider  an  explosive  source  centered  at  the  origin  and  assume  cotnpleie  .spherical  symineuy  so 
that  all  displacements  will  be  in  the  radial  direction  and  only  coinpre.ssioiial  waves  will  be  generated. 
Adopting  a  spherical  coordinate  system  and  assuming  a  harmonic  time  dependence  of  tiic  lonn  tlie 
displacement  can  be  written  as 

u  =  u(r,(o)r  (4.1) 


where  f  is  the  unit  vector  in  the  radial  direction.  In  a  region  where  the  equtitions  ol  linctir  elasticity  ;tre 
obeyed  the  displacement  must  satisfy 

(X.+2|i)V^u  -  VxVu  +  pto'u  =  0  (4.2) 

where  p  is  the  density  and  X  ,  p  are  the  Lame  elastic  constants.  In  a  regioit  where  these  elastic  pttnun- 
eters  are  independent  of  the  spatial  coordinates,  a  general  solution  to  the  tihovc  cquaiion  is 

ufr.o))  =  +  a<^^/t{^Uipr)r  (4.}) 

Here  and  are  spherical  Hankel  functions  of  the  first  and  .second  kind,  respectively,  and 


(4,4) 


where 


It  is  easily  demonstrated  that  (he  function  A/”  repre-sents  inward  traveling  coinpressiimal  wa\es,  while 
represents  outward  traveling  waves,  'fhe  terms  rr*'*  and  «*•’  :wc  arbitrary  consianis  ;ii  this  point. 

The  radial  traction  on  any  surface  of  constant  radius  has  the  lonn 


XV  ur  + 


u  +  pfxVxu  =  x„(r.(o)f 


(4,6) 


Corre,sponding  to  the  solution  for  u  given  above,  this  stress  is 


trr(r,ci))  =  a 


(1) 


(X+2p)kpAo'’  (k,,  r) - ^A  (tpr ) 


+  a 


<2) 


(X+2p )<,•,,  A/,-’  (kp  r )  -  ^ A  <1).  r ) 


t4.7) 
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The  tail  ,  .iial  stresses  arc 


Xee(r.«B)  =  t^Cr.co)  =  a 


+  a 


(2) 


'Kkphi^\k.T)  +  ^h{^Hk,r) 


(4.8) 


Now  consider  the  region  around  the  source  to  be  divided  into  a  number  of  spherical  shells,  with 
the  material  properties  constant  within  each  shell  but  possibly  changing  at  the  boundiuies  between 
^lls.  Shell  n  will  be  bounded  by  radii  r,  and  r^+i  and  will  have  material  properties  p„.  X„,  p„.  We 
also  have  k„  =  m/Vp,, .  Using  the  results  developed  above,  the  dlsplaceinem  ;uid  radial  sUcss  wiUiin  the 
n-th  shell  can  be  written  in  matrix  form  as 


u 

tr 


=  A, 


(4.<J) 


where  the  shell  matrix  A  is  defined  as 


A,(r,a) 


/I  I-' 


Uv) 


4^, 


4m„ 


[(K-^-2\in)k^h^^Hk,r)  -  -r^h{^Hk„r)]  [(X,+2iiJit„ {f„r)  - 


(4.10) 


It  is  easily  shown  that  the  inverse  of  this  matrix  exists  and  is  given  by 


A;‘(r,co)  = 


ik,r^ 


2(X,+2fi„) 


l(X„+2p„)I-,h,«>(I:„r)-  ^/ip(I„r)]  -hfHk„r) 
-[(K^2\i„)k„h(,'\k„r)  -  ^l,i'>{k,.r)]  fil‘>(k.,r) 


(4.11) 


The  boundary  conditions  for  this  problem  are  that  u  and  f,  should  be  everywherc  continuous. 
This  means  that  at  the  radius  r„  where  shells  n  and  n-l  are  in  conuict  die  following  ciiuation  must  hold. 

,<1)1 


=  A,_j(r„,co) 


=  A„(r„,©) 


,(2) 


(4.12) 


It  follows  that 

,(i)' 


fl(2) 

L“"  J 


=  A;‘(r,,fi>)  A„.i(r„.(a) 


„(i) 

“o-l 

«(2) 


(4.1,^) 


It  is  convenient  to  define  an  interface  matrix  as 


B„(®)  s  A;'(r„,to)  A„-,(r„,(i)) 
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(4.14) 


Letting  z„  -  k„r„  and  :„_i  =  the  individual  elements  of  B„  are 


Bu  = 


= 


z„U  z-  z^ZnU 


2(X,+2h„) 


^  ^  *  z«-l  \  ^  ^  A  ■  /  \  Z^Ki  Jn-l) 

(X,+2n, )— ^5 - a„.,+2M,.,)— y-+4/  (^,  - — - 

Z--1  Z.  Z,Z„*I 


(4,15b) 


(4.150 

2(^+2li,)  z„-.i  z/  z„h„.i 


Bu  =  TrirH.^M  r  •*'2lt.. )-  -(X,.i+2M^.,)-^+4i' -n„_i)  "  yV"'‘^  (4.15d) 

/  Z#i-t  Zfl  J 

This  process  of  matching  boundary  conditions  at  an  interface  can  be  repeated  at  succe.ssive  boundtuies, 
and  thus  the  solution  can  be  propagated  through  an  arbitrary  number  of  shells.  A.ssuming  the  process  to 
have  started  in  shell  1,  the  solution  in  shell  n  is 


,,U)  =  B,(o))  Bold))  ,2) 


The  next  step  is  to  assume  that  shell  0  consists  of  a  gas  having  pressure  P  and  adiabatic  gas  con¬ 
stant  Y-  Then,  allowing  approximately  for  the  change  in  pre.ssure  caused  by  the  exp.ansion  of  tlie  shell, 
the  displacement-stres.s  vector  at  the  radius  will  be 

u  u 

r,j  =  [-P(l-3YU/ri)  =  |^^(2)J  (•^■17) 

Two  independent  trial  solutions  can  be  ohitiined  for  Uiis  system.  First  a.s,sumc  that  « {*'  =  0  and  theti  it 
is  possible  to  solve  for 

(o  -P  . 


The  second  solution  is  obtained  by  assuming  that  «{'*  =  0  ;uid  then 
..  t2)  -P 


(4.18a) 


ap  = 


(X,+2^,)*,/2, r,,  -  — (n,  +  ^)/tp(f-,r,) 
Ti  4 


(4.18b) 


The  general  solution  is  obtained  by  introducing  the  p:ir;uncter  E,  and  tlicn  tltc  sum  of  E  times  the  first 
trial  solution  :uid  (1-e)  limes  the  second  trial  solution  will  .sati.sly  the  hound;u7  conditions  at  r,. 


1.3 


Tbe  parameter  e  is  determined  by  the  radiation  condition.  Assume  iliat  shell  N  extends  to  inliniiy 
so  that  there  are  only  outward  propagating  wave.s  in  this  region,  and  thus  Uji/'  ’  =  0.  I'his  inemis  tliat 

,  EfliiV  +  a-e)a/^V  =  0  *  (4,1‘J) 


where  Ujvi  and  are  the  two  solutirxis  that  are  obtained  in  .shell  N  when  the  two  dilicreni  iruil  .solu 
tions  from  rj  are  propagated  independently  to  r//.  It  follows  that 


E  = 


(4.20) 


Having  determined  e,  the  two  trial  solutions  can  be  summed  with  the  proper  weighting  laeiors  at  tmy 
radius  to  yield  the  total  solution. 

The  case  of  an  explosion  in  a  medium  where  a  hydrostatic  pressure  is  present  is  obuiined  rn)m 
the  above  solution  by  letting  the  pressure  on  the  cavity  wall  be  P  -  P/ ,  where  P/  is  die  lithostatie  pres¬ 
sure  at  the  shot  point.  Then  all  of  the  pressures  and  suesses  that  are  ctilcuhited  will  be  departures  Iroin 
the  initial  hydrostatic  case.  With  this  approach  the  initial  stre.ss  t„  =  -P^I  is  homogeneous  so  V  t„  =  0 
and  there  is  no  change  is  the  basic  equations.  However,  if  the  initial  pressure  is  not  hoinogcneous,  such 
as  when  the  depth  effect  of  pressure  is  included,  the  .solution  loses  its  spherical  symmetry  and  becomes 
more  complicated. 

The  possibility  that  elastic  waves  will  be  attenuated  a.s  they  propagate  is  itieludcd  by  introducing 
anelasticity  into  the  problem.  This  is  achieved  by  u.sing  the  correspondence  principle  to  generali/.e  die 
previous  results.  This  consists  of  introducing  the  complex  elastic  coiisttuiis 

X  +  2p  -»  (X  +  2p)(l  +  iQ-'}  (4.21a) 

and 


p  n(l  +  iQr‘)  (4.21b) 

The  wave  number  ip  will  now  be  complex  also  and  can  be  written 

ip  lfcp(l-i8)  (4.22) 

The  requirement  that  the  elastodynamic  equations  be  satislied  leads  to 

•/i 

(4.2".a) 

and 

5  =  - - - nr 

Qpd  +  [!+(?, r-1^) 


CO 


1 


1  +  [i+<2; 


x+2p  1+a 


«2 
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A  similar  set  of  equations  can  be  derived  for  the  wavenumber  A',.  Note  that  in  general  and  Q,  w;;; 
be  functions  of  frequency. 


With  attenuation  present,  the  vekx;ities  will  in  general  be  dispersive.  Assuming  that  tind  (J, 
are  not  functions  of  frequency  in  the  band  of  interest,  cau.saliiy  is  mtiintaihed  by  requiring  that 


Vp(f)  =  V,(f,) 

Vsif)  =  y.fjr) 


1  +  -^ln[-f  ] 

llQp  fr 


1  + 

^Qs  fr 


{4.24a) 


(4,24b) 


where  is  some  reference  frequency. 


5.  Method  of  Calculation 

The  basic  steps  of  the  computational  scheme  are  outlined  in  Figure  3.  The  lirsi  step  is  to  inititil- 
ize  the  model  of  the  material  properties  surrounding  the  explo.sion.  This  consists  of  dividing  tlie  region 
into  an  appropriate  number  of  spherictU  shells  and  specifying  the  inititil  material  properties  for  ihe.se 
shells.  The  material  properties  tue  assumed  to  he  con.stant  within  each  shell.  The  piiramciers  associ¬ 
ated  with  the  equivalent  elastic  treatment,  such  as  the  ultimate  density,  the  reference  shetir  strain,  ;uid 
the  maximum  damping  must  al.so  be  specilied  at  this  time.  The  outside  radius  of  the  last  shell  should 
be  large  enough  so  that  the  material  properties  at  tliis  and  greater  disttuices  can  be  stilcly  tissumed  to  be 
linear  visco-elastic. 

The  second  step  is  to  initialize  the  properties  of  the  explosive  source.  The  size  of  the  explosion 
is  .specified  by  giving  its  chemical  energy  £„.  Tlie  initial  radius  of  tlie  ctivity  r„  is  al.so  specilied. 
which  leads  to  the  initial  volume 


The  explo.sion  is  assumed  to  Uike  place  instantancou.sly  .so  that  at  t  =  0  the  initial  pressure  within  die 
cavity  increases  by  an  tunount 


V.. 


(.s.:) 


where  y  is  Uic  gas  consltmt  for  the  explosive  gas,sc.s.  This  inititil  pressure  is  the  radial  stress  at  r  = 
which  is  one  of  tlie  boundtiry  conditions  for  the  elastic  wave  .solution. 
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Figure  3.  Flow  chart  for  the  calculation  ol  the  outward  propagating  waves  fn^m  an  explosive  source 
using  the  equivalent  elastic  methrxl. 
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The  third  step  is  to  obtain  the  elastic  wave  srdution  for  the  entire  region  surrounding  the  source 
using  the  propagator  solutions  of  section  4.  Given  this  solution,  the  stres.ses  and  strains  at  ;dl  distmices 
from  the  source  can  be  “calculated.  In  generiU  the  su-es.ses  as.stK:iated  with  this  solution  will  be 
sufficiently  large  so  that  the  initial  material  properties  are  not  appropriate.  Thus  the  material  properties 
are  adjusted  on  the  basis  of  the  calculated  stresses  and  the  equivalent  elastic  u-eatment  outlined  in  sec¬ 
tion  2.  For  the  purposes  of  tliis  adjustment  the  pressure  is 


P  =  - 


X,,  +x^  +  x^ 


and  the  shear  stfess  is 


Having  adjusted  the  material  properties  in  all  of  the  shells,  the  elastic  wave  solution  c;in  he  obtained  for 
the  adjusted  mcxlel  and  the  entire  pr(x:ess  repeated.  This  iterative  process  of  solving  for  the  stresses  ;uid 
adjusting  the  material  properties  is  repeated  until  a  solution  is  obtained  in  which  tfic  stresses  are  compa¬ 
tible  with  the  material  properties  at  all  distances.  Generally,  only  three  or  four  iterations  ;ue  required 

After  convergence  of  the  prix.’ess  of  adjusting  tire  material  properties,  the  source  radius  is  also 
adjusted  to  allow  for  the  inelastic  growth  of  the  cavity.  The  cavity  radius  r„  is  increased  by  succes¬ 
sively  eliminating  shells  until  tlie  static  pressure  is  less  th;ui  some  factor  times  the  litluvstatic  pressure 

P(r„ .!-><«)  <  n  P,  {5.5) 

Here  the  lilhostatic  pressure  P/  is  just  ditit  due  to  the  overburden  and  is  given  by 

Pi  =  Pr,  g 

where  g  is  the  acceleration  of  gravity  ;uid  is  the  .source  dcpili.  Recall  dial  llie  calculated  prcssuie  is 
in  excess  of  the  lilhostatic  pressure,  .so  diis  condition  actually  corresponds  to  the  static  cavity  pressure 
being  l-t-T|  times  the  lilhostatic  pre.ssure.  For  the  calculations  in  luff  die  factor  ri  was  icssumcd  to  have 
a  value  of  0.5. 

After  die  adjustment  of  the  cavity  radius,  a  tinaJ  iteration  is  performed  on  the  adjusimenl  ol  die 
material  properties  to  insure  that  the  sire.s.ses  and  material  properties  at  ;dl  distances  from  die  source  arc 
compatible.  At  die  completion  of  this  step  die  final  .solution  is  available. 


6.  Expt:rim«ntal  Check  uf  the  Method 

This  method  of  calculating  fields  of  an  explosive  source  is  illustrated  by  some  initial  calculations 
in  Figure  4.  The  upper  pimel  shows  the  first  10  msec  of  measured  particle  veftK'iiy  within  a  lew  meters 
of  a  buried  chemical  explosion,  while  die  lower  panel  shows  the  results  of  simulating  these  measure¬ 
ments  with  the  code  that  employs  equivalent  elastic  material  properties.  The  data  were  obuiined  during 
the  OSSY2  experiments  of  1991. 

The  OSSY2  (On  Site  Seismic  Yield)  experiments  were  perlonned  in  1991  in  Yucca  Valley  at  die 
Nevada  Test  Site.  They  were  a  cooperative  effort  between  scientists  at  Lawrence  Livermore  National 
Laboratory,  Lawrence  Berkeley  Laboratory,  the  Seismograpliic  Suition  at  UC  Berkeley,  and  Southern 
Methodist  Llniversity.  The  dam  shown  in  Figure  4  were  obtained  from  a  chemical  explo.sion  dial  was 
detonated  at  a  depth  of  534  meters  in  hole  L)E4av.  The  .source  size  was  100  pounds  ot  ('4  explosive 
The  source  medium  was  a  bedded  tulf  immediately  above  the  Grouse  Canyon  member.  The  source 
properties  which  were  used  in  the  calculation  are 

p,  =  2.10  g«i/cnr’  ,  Vp,  =  2.^5  km  I  sec  ,  V^„  =  1.45  Ani/.vcc 

As  mentioned  in  section  2,  the  equivalent  cktsiic  parameters  u.sed  in  the  calculations  were 
Pu  =  2.58gwi/cvn  \  =  i0~*,  and  dintix  =  0.30.  These  estimates  ot  material  properties  were  based  on 

the  re.sults  of  a  VSP  survey  (Leonitrd  et  al,  1992)  ;uid  laboratory  measurements  (Fanerick.  1966). 

Starting  with  the  material  properties  listed  above  as  initial  vidoes,  the  el.astic  wave  solutions  were 
obtained  out  to  a  distance  of  40  meters  from  the  source  point.  Using  the  stresses  diat  resulted  Ifoin  diis 
solution,  the  material  properties  were  adjusted  according  to  the  mctliod  described  in  section  2.  I'hcii  a 
new  wave  solution  was  obtained  using  the  modified  material  properties.  This  process  was  repeated  until 
there  was  essentially  no  change  in  material  properties  between  iienitions. 

The  equivident  clastic  modification  uf  the  material  properties  was  most  pronounced  in  a  region 
that  extended  out  to  about  2  meters  from  llie  shot  point.  In  this  region  the  density  and  coinpressiunal 
velocity  were  increased  by  up  to  5%  while  the  shciu  velocity  ;uid  shear  quiUily  factor  were  decreased 
by  over  90%.  At  liager  distances  there  was  little  modilication  in  die  density  and  coinpre.ssional  velo¬ 
city,  although  there  were  signilicant  but  decreasing  elfecls  upon  die  slietir  velocity  and  siicar  qtialiiy 
factor  out  to  disiiuiccs  of  about  10  meters. 
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Figure  4.  Recorded  and  calculated  free  field  measurements  of  velocities  from  a  chemical  explosion. 
The  explosion  consisted  of  100  pounds  of  C4  explosive  which  was  detonated  at  a  depth  of  5M  meters 
in  Yucca  Valley.  The  upper  panel  shows  velocities  that  were  recorded  by  a  linear  array  of  instruments 
directly  above  the  explosion.  The  number  on  the  left  of  each  trace  is  the  radial  distmicc  of  the  instru¬ 


ment  from  the  center  of  the  explosion  in  meters  and  the  number  on  the  right  is  the  maximum  vehKity 
in  cm/sec.  The  traces  have  been  multiplied  by  distance  Irom  tlie  explosion  for  the  purpo.scs  of  pUttting, 
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The  next  step  was  to  calculate  a  new  en'ecii%'e  cavity  radius  for  the  source.  The  criterion  used 
was  that  the  static  cavity  pressure  should  equal  1.5  times  the  lithosiatic  pressure.  This  resulted  in  die 
cavity  radius  increasing  from  9  cm  to  20  cm.  Another  iteration  of  the  equiviilent  eltisiic  adjustment  of 
material  properties  was  then  perfonned,  with  the  main  result  being  a  further  reduction  in  the  shetir  velo¬ 
city  in  the  immediate  vicinity  of  the  new  cavity  radius.  Finally,  using  the  effective  material  properties 
and  the  effective  cavity  radius,  a  fmtil  wave  solution  was  obtained  to  yield  the  results  shown  in  the 
lower  panel  of  Figure  4. 

The  agreement  between  the  observed  and  calculated  velocities  in  Figure  4  indicates  .some  promi.se 
for  this  type  of  an  approach.  The  recordings  at  the  two  closest  distances  of  1.5  and  7.2  meters  ;ae 
somewhat  suspect,  so  the  comparison  here  may  not  be  valid.  The  acceleration  at  the  first  gage  was  m 
excess  of  8000  g  and  there  is  a  suggestion  on  the  acceleration  record  that  the  gage  may  have  broken 
kxxse  and  gone  into  ballistic  motion.  The  acceleration  at  the  second  gage  was  greater  th;ui  251X)  g  ;uid 
tlie  maximum  acceleration  iKcurs  during  the  second  pul.se  that  is  delayed  by  about  .1  msec,  which  is 
difficult  to  explain  in  terms  of  an  outgoing  pre.ssure  wave.  Note  that  both  of  these  gages  were  in  the 
range  of  su-ongly  nonlinear  material  behavior  in  the  sense  that  the  equivalent  clastic  u-eauneni  resuiied 
in  significant  mtxlifications.  At  the  three  outer  gages  the  observed  records  .seem  more  reasontihle.  indi¬ 
cating  an  outward  propagating  wave  that  chtuigcs  slowly  with  distance.  The  simulated  records  agree 
quite  well  in  amplitude  ;uid  period  in  this  nuige,  although  die  asymtneU'y  in  die  waveforms  in  somewhat 
different  for  the  observed  and  etdculaied  results.  One  ptissihle  expkuiatioii  for  this  difference  is  (h;ti  die 
dispersion  asstKiated  with  the  anelastie  properties  of  die  medium  has  not  iieen  properly  modeled  in  die 
.simulations.  Tlii.s,  tilting  with  intuiy  oilier  aspects  of  the  simultition  calcuhiiions.  need  considerably 
more  investigation. 

7.  Conclusions 

The  approach  dc.scribed  here  is  not  a  substitute  for  complete  equatioii-of-st;itc  codes,  hut  it 
appears  to  provide  ;ui  effective  mediod  of  investigating  .some  of  die  basic  elements  tif  wave  propagation 
iietu'  tin  explosive  source,  'llic  computer  codes  tire  quite  efliciciu,  so  that  it  is  ptissihlc  to  pcifonn 
extensive  siinulation  studies  in  which  the  role  of  various  ptutuncters  is  investigated.  Thus  die  tnelhod  is 
quite  u.seful  in  helping  to  isolate  die  model  partunclers  which  arc  controlling  sonic  ptirticular  fetiiuie  of 
the  observational  dtilti. 

Considerable  more  testing  must  be  perfonned  in  order  to  establish  Unit  this  lucthod  piodiices  ivli- 
ablc  qutintilative  estimatc.s  of  the  seismic  wtivcforms  which  arc  traiiMiiiiied  to  the  linetir  elasite  region. 
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So  far  most  of  the  comparisons  with  observatiomil  data  have  involved  chemical  explosions,  because 
.such  data  were  conveniently  available,  but  such  compmsons  must  be  extended  to  nuclciir  explosions.  It 
would  also  be  helpful  to  extend  the  comptirisons  to  other  types  of  source  media. 

One  characteristic  of  the  equivalent  elastic  method  is  that  it  requires  a  rehuionship  between 
material  properties  and  the  level  of  strain.  Further  deveUtpment  and  use  of  this  method  will  require  the 
assembly  of  a  data  ba.se  of  em-'irical  relationships  for  equivalent  elastic  treatments  that  covers  till  of  the 
different  types  of  materials  which  might  serve  as  source  media  for  explosions.  The  liteniture  conuiins 
considerable  data  of  this  type  for  various  soils,  but  comptuable  dtita  pertaining  to  h;ad  rocks  htts  not  yet 
been  ItKated.  However,  a  complete  search  of  the  literature  with  respect  to  this  type  of  information  has 
not  been  completed. 

While  the  feasibility  of  the  method  investigated  in  this  .study  has  iu)t  yet  been  tinnly  establi.shcd. 
the  exercise  has  been  successful  in  demonstrating  the  need  for  approaches  of  this  type  which  produce 
simple  analytic  results.  One  of  the  advtuitages  of  such  tin  approach  is  thtit  it  tUlows  vtirious  lealures  of 
the  observational  data  to  be  explored  with  re.spect  to  ilieir  cau.se  tutd  uniqueness.  The  reltitive  role  of 
the  material  inside  and  outside  the  eltLstic  radius  in  determining  the  frequency  content  of  die  nidiaied 
.seismic  energy  is  one  such  area  of  interest.  Another  is  die  relationship  between  tlie  citnier  frcqucticy  of 
tlie  spectrum  and  the  source  radius.  For  in.sttince.  the  itiitial  calcuhitions  suggest  thtu  tltc  clfective 
.source  radius  may  be  a  rather  vaguely  defined  quantity.  It  is  also  po.ssible  to  explore  Uie  role  of  die 
energy  density  within  the  cavity,  which  relates  to  die  differences  lieiween  cheinictil  and  nuclettr  explo¬ 
sions  and  methods  of  decoupling  explosions  in  ;u)  cnhu'gcd  cavity. 
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